GIS 713 | Preliminary Analyses
G. Millar
Setup
library(sp)
library(raster)
library(R.matlab)
library(plyr)
library(data.table)
library(maptools)
library(rgdal)
library(spatstat) Data Proccessing
setwd("~/Desktop/Fall_2018_Classes/GIS_713/Final_Project")
CHIPS_df <- read.table("chips.csv", header = TRUE, row.names=NULL, sep=",")
CHIPS_df <- CHIPS_df[!CHIPS_df$lat < 500000.00, ]
CHIPS_df <- CHIPS_df[!CHIPS_df$long < 50000.00, ]
# Correcting odd lat coords
CHIPS_df$lat <- CHIPS_df$lat / 10000
CHIPS_df$long <- CHIPS_df$long / 10000
# And a function to shift vectors conviniently:
shift.vec <- function (vec, shift) {
if(length(vec) <= abs(shift)) {
rep(NA ,length(vec))
}else{
if (shift >= 0) {
c(rep(NA, shift), vec[1:(length(vec)-shift)]) }
else {
c(vec[(abs(shift)+1):length(vec)], rep(NA, abs(shift))) } }
}
# Calculating distances between successive positions and the respective speed in this segment.
# Shifting vectors for lat and lon so that each row also contains the next position.
CHIPS_df$lat.p1 <- shift.vec(CHIPS_df$lat, -1)
CHIPS_df$lon.p1 <- shift.vec(CHIPS_df$long, -1)
# Calculating distances (in metres) using the function pointDistance from the ‘raster’ package.
CHIPS_df$dist.to.prev <- apply(CHIPS_df, 1, FUN = function (row) {
pointDistance(c(as.numeric(row["lat.p1"]), as.numeric(row["long.p1"])),
c(as.numeric(row["lat"]), as.numeric(row["long"])),
lonlat = T)
})
# Transforming the column ‘time’ so that R knows how to interpret it.
CHIPS_df$time_new <- strptime(CHIPS_df$initial_time_stamp_mat,
format="%m/%d/%Y %H:%M")
# Shift the time vector, too.
CHIPS_df$time.p1 <- shift.vec(CHIPS_df$time_new, -1)
# Calculating number of seconds between two positions.
CHIPS_df$time.diff.to.prev <- as.numeric(difftime(CHIPS_df$time.p1,
CHIPS_df$time_new))
# Calculating metres per seconds, kilometres per hour, and two LOWESS smoothers to get rid of some noise.
CHIPS_df$speed.m.per.sec <- CHIPS_df$dist.to.prev / CHIPS_df$time.diff.to.prev
CHIPS_df$speed.km.per.h <- CHIPS_df$speed.m.per.sec * 3.6
CHIPS_df$speed.km.per.h <- ifelse(is.na(CHIPS_df$speed.km.per.h), 0,
CHIPS_df$speed.km.per.h)
CHIPS_df$lowess.speed <- lowess(CHIPS_df$speed2, f = 0.2)$y
CHIPS_df$lowess.alt <- lowess(CHIPS_df$altitude, f = 0.2)$y
CHIPS_df$lowess.conduct <- lowess(CHIPS_df$conductance_z, f = 0.2)$y
setnames(CHIPS_df, "long", "lon")
pt1 <- CHIPS_df[CHIPS_df$participant == 1, ]
pt2 <- CHIPS_df[CHIPS_df$participant == 2, ]
pt3 <- CHIPS_df[CHIPS_df$participant == 3, ]
pt4 <- CHIPS_df[CHIPS_df$participant == 4, ]
pt5 <- CHIPS_df[CHIPS_df$participant == 5, ]
pt6 <- CHIPS_df[CHIPS_df$participant == 6, ]
pt7 <- CHIPS_df[CHIPS_df$participant == 7, ]
pt8 <- CHIPS_df[CHIPS_df$participant == 8, ]
pt9 <- CHIPS_df[CHIPS_df$participant == 9, ]
pt10 <- CHIPS_df[CHIPS_df$participant == 10, ]
pt11 <- CHIPS_df[CHIPS_df$participant == 11, ]Exploratory Data Visualizations
GPS
# Now, let’s plot all the stuff!
# Plot elevations and smoother
layout(matrix(1:3, nrow=3))
plot(CHIPS_df$altitude, type = "l", bty = "n", xaxt = "n", lwd= 3,
ylab = "Elevation",
xlab = "", col = "grey60")
lines(CHIPS_df$lowess.alt, col = "green", lwd = 3)
abline(h = mean(CHIPS_df$altitude), lty = 2, lwd = 3, col = "green")
legend(x="bottomright", legend = c("GPS elevation", "LOWESS elevation",
"Mean elevation"),
col = c("grey60", "green", "green"), lwd = c(2,4,2), lty = c(1,2,2),
bty = "n")
# Plot speeds and smoother
plot(CHIPS_df$speed2, type = "l", bty = "n", lwd= 3, xaxt = "n",
ylab = "Speed (km/h)", xlab = "", col = "grey60")
lines(CHIPS_df$lowess.speed, col = "red", lwd = 3)
abline(h = mean(CHIPS_df$speed2), lty = 2, lwd = 3, col = "red")
legend(x="topright", legend = c("GPS speed", "LOWESS speed",
"Mean speed"),
col = c("grey60", "red", "red"), lwd = c(2,4,2), lty = c(1,2,2),
bty = "n")
# Plot conductnace and smoother
plot(CHIPS_df$conductance_z, type = "l", bty = "n", lwd= 3, xaxt = "n",
ylab = "Skin Conductance", xlab = "", col = "grey60")
lines(CHIPS_df$lowess.conduct, col = "blue", lwd = 3)
abline(h = mean(CHIPS_df$conductance_z), lty = 2, lwd = 3, col = "blue")
legend(x="topright",
legend = c("Conductance", "LOWESS conductance", "Mean conductance"),
col = c("grey60", "blue", "blue"), lwd = c(2,4,2), lty = c(1,2,2),
bty = "n")
par(mar=c(5, 4, 4, 2) + 0.1)Skin Conductance
Plotting the elevation and timestamp of each waypoint using ggplot, allows the visualization of the cycling route between Tilburg and Waalwijk. As a preparatory step, the ymd_hms() function from the lubridate library is used to convert the string representating the timestamp to a proper R time-object. To not confuse ggplot, the SpatialPointsDataFrame-object is not passed directly, but converted to a regular dataframe with as.data.frame():
if(!require(lubridate)) install_github("rstudio/lubridate")
if(!require(ggplot2)) install_github("rstudio/ggplot2")
if(!require(gridExtra)) install_github("rstudio/gridExtra")
# plot of time and elevation, colored by skin conductance
time_ele_conduct_plot <- ggplot(as.data.frame(CHIPS_df), # convert to regular dataframe
aes(x=time, y=altitude, color = conductance_z)) +
scale_color_gradient(low="blue", high="red") +
geom_point(alpha = 0.01, size = 2) +
labs(x='Cycling time', y='Elevations (meters)')
# plot of time and speed, colored by skin conductance
time_speed_conduct_plot <- ggplot(as.data.frame(CHIPS_df), # convert to regular dataframe
aes(x=time, y=speed2, color = conductance_z)) +
scale_color_gradient(low="blue", high="red") +
geom_point(alpha = 0.08, size = 2) +
labs(x='Cycling time', y='Speed (km/h)')
grid.arrange(time_ele_conduct_plot, time_speed_conduct_plot, nrow=2)
Spatial Data Processing
shp_dsn <- "~/Desktop/Fall_2018_Classes/GIS_713/Final_Project/NL006L3_TILBURG/Shapefiles"
landcover <- readOGR(path.expand(shp_dsn), 'NL006L3_TILBURG_UA2012')## OGR data source with driver: ESRI Shapefile
## Source: "/Users/garrettmillar/Desktop/Fall_2018_Classes/GIS_713/Final_Project/NL006L3_TILBURG/Shapefiles", layer: "NL006L3_TILBURG_UA2012"
## with 12086 features
## It has 11 fields# Projection
landcover <- spTransform(landcover, CRS("+proj=longlat +datum=WGS84"))
# Conversion into SpatialPoints
coordinates(CHIPS_df) <- ~ lon + lat
coordinates(pt1) <- ~ lon + lat
coordinates(pt2) <- ~ lon + lat
coordinates(pt3) <- ~ lon + lat
coordinates(pt4) <- ~ lon + lat
coordinates(pt5) <- ~ lon + lat
coordinates(pt6) <- ~ lon + lat
coordinates(pt7) <- ~ lon + lat
coordinates(pt8) <- ~ lon + lat
coordinates(pt9) <- ~ lon + lat
coordinates(pt10) <- ~ lon + lat
coordinates(pt11) <- ~ lon + lat
# Setting default projection
proj4string(CHIPS_df) <- CRS('+proj=longlat +datum=WGS84')
proj4string(pt1) <- CRS('+proj=longlat +datum=WGS84')
proj4string(pt2) <- CRS('+proj=longlat +datum=WGS84')
proj4string(pt3) <- CRS('+proj=longlat +datum=WGS84')
proj4string(pt4) <- CRS('+proj=longlat +datum=WGS84')
proj4string(pt5) <- CRS('+proj=longlat +datum=WGS84')
proj4string(pt6) <- CRS('+proj=longlat +datum=WGS84')
proj4string(pt7) <- CRS('+proj=longlat +datum=WGS84')
proj4string(pt8) <- CRS('+proj=longlat +datum=WGS84')
proj4string(pt9) <- CRS('+proj=longlat +datum=WGS84')
proj4string(pt10) <- CRS('+proj=longlat +datum=WGS84')
proj4string(pt11) <- CRS('+proj=longlat +datum=WGS84')Web-mapping
library(leaflet)
require(pals)
conduct.pal <- colorNumeric (c("dodgerblue4", "slategray2", "red3"),
pt1$conductance_z)
m <- leaflet() %>%
# Add tiles
addProviderTiles("Esri.WorldTopoMap", group = "Topographical") %>%
addProviderTiles("OpenStreetMap.Mapnik", group = "Road map") %>%
addProviderTiles("Esri.WorldImagery", group = "Satellite") %>%
addCircles (data=pt1, group='Participant 1', stroke = T, radius = 80,
fillOpacity = 0.2, fillColor = conduct.pal(pt1$conductance_z),
opacity = 0.2, color = conduct.pal(pt1$conductance_z)) %>%
addCircles (data=pt2, group='Participant 2', stroke = T, radius = 80,
fillOpacity = 0.2, fillColor = conduct.pal(pt2$conductance_z),
opacity = 0.2, color = conduct.pal(pt2$conductance_z)) %>%
addCircles (data=pt3, group='Participant 3', stroke = T, radius = 80,
fillOpacity = 0.2, fillColor = conduct.pal(pt3$conductance_z),
opacity = 0.2, color = conduct.pal(pt3$conductance_z)) %>%
addCircles (data=pt4, group='Participant 4', stroke = T, radius = 80,
fillOpacity = 0.2, fillColor = conduct.pal(pt4$conductance_z),
opacity = 0.2, color = conduct.pal(pt4$conductance_z)) %>%
addCircles (data=pt5, group='Participant 5', stroke = T, radius = 80,
fillOpacity = 0.2, fillColor = conduct.pal(pt5$conductance_z),
opacity = 0.2, color = conduct.pal(pt5$conductance_z)) %>%
addCircles (data=pt6, group='Participant 6', stroke = T, radius = 80,
fillOpacity = 0.2, fillColor = conduct.pal(pt6$conductance_z),
opacity = 0.2, color = conduct.pal(pt6$conductance_z)) %>%
addCircles (data=pt7, group='Participant 7', stroke = T, radius = 80,
fillOpacity = 0.2, fillColor = conduct.pal(pt7$conductance_z),
opacity = 0.2, color = conduct.pal(pt7$conductance_z)) %>%
addCircles (data=pt8, group='Participant 8', stroke = T, radius = 80,
fillOpacity = 0.2, fillColor = conduct.pal(pt8$conductance_z),
opacity = 0.2, color = conduct.pal(pt8$conductance_z)) %>%
addCircles (data=pt9, group='Participant 9', stroke = T, radius = 80,
fillOpacity = 0.2, fillColor = conduct.pal(pt9$conductance_z),
opacity = 0.2, color = conduct.pal(pt9$conductance_z)) %>%
addCircles (data=pt10, group='Participant 10', stroke = T, radius = 80,
fillOpacity = 0.2, fillColor = conduct.pal(pt10$conductance_z),
opacity = 0.2, color = conduct.pal(pt10$conductance_z)) %>%
addCircles (data=pt11, group='Participant 11', stroke = T, radius = 80,
fillOpacity = 0.2, fillColor = conduct.pal(pt11$conductance_z),
opacity = 0.2, color = conduct.pal(pt11$conductance_z)) %>%
# Layers control
addLayersControl(position = 'bottomright',
baseGroups = c("Topographical", "Road map", "Satellite"),
overlayGroups = c("Participant 1", "Participant 2",
"Participant 3", "Participant 4",
"Participant 5", "Participant 6",
"Participant 7", "Participant 8",
"Participant 9", "Participant 10",
"Participant 11"),
options = layersControlOptions(collapsed = FALSE)) %>%
hideGroup(c("Participant 2", "Participant 3", "Participant 4", "Participant 5",
"Participant 6", "Participant 7", "Participant 8", "Participant 9",
"Participant 10", "Participant 11")) %>%
addLegend(values = pt1$conductance_z, pal = conduct.pal,
opacity = 1, title = "Skin Conductivity", position = "bottomleft")
mStudy Area
urban <- c("Discontinuous dense urban fabric (S.L. : 50% - 80%)",
"Industrial, commercial, public, military and private units",
"Discontinuous low density urban fabric (S.L. : 10% - 30%)",
"Isolated structures",
"Continuous urban fabric (S.L. : > 80%)",
"Discontinuous very low density urban fabric (S.L. : < 10%)",
"Discontinuous medium density urban fabric (S.L. : 30% - 50%)",
"Construction sites",
"Mineral extraction and dump sites")
natural <- c("Arable land (annual crops)",
"Pastures",
"Forests",
"Land without current use",
"Green urban areas",
"Permanent crops (vineyards, fruit trees, olive groves)",
"Sports and leisure facilities",
"Herbaceous vegetation associations (natural grassland, moors...)",
"Open spaces with little or no vegetation (beaches, dunes, bare rocks, glaciers)",
"Wetlands")
roads <- c("Other roads and associated land",
"Railways and associated land",
"Fast transit roads and associated land")
name = landcover$ITEM2012
landcover$group <- with(landcover, ifelse(name %in% urban, "urban",
ifelse(name %in% natural,
"natural", "roads")))
col = terrain.colors(3)
plot(landcover, col = col, col.regions = landcover$group,
edge.col = "transparent", axes = F,
colorkey = list(space = "bottom", height = 0.5,
width = 0.7),
main = "Study Area", sub = "Types of Landcover",
par.settings = list(axis.line = list(col = 'transparent')))
legend("bottomleft", title = "Landcover", fill = col,
legend = c("Urban", "Natural", "Roads"))
Spatial Analysis
Mapping Cycling Routes
# clipping landcover polygon to cycling route
neth_clipped <- as(crop(landcover, CHIPS_df), "SpatialPolygonsDataFrame")
col = terrain.colors(3)
conduct.pal <- colorNumeric (c("dodgerblue4", "slategray2", "red3"),
pt1$conductance_z)
layout(matrix(1:6, nrow=1))
plot(neth_clipped, col = col, col.regions = neth_clipped$group,
edge.col = "transparent", axes = F,
colorkey = list(space = "bottom", height = 0.5, width = 0.7),
main = "",
par.settings = list(axis.line = list(col = 'transparent')))
points(pt1, col = conduct.pal(pt1$conductance_z), cex = 1.5, pch = 20 )
plot(neth_clipped, col = col, col.regions = neth_clipped$group,
edge.col = "transparent", axes = F,
colorkey = list(space = "bottom", height = 0.5, width = 0.7),
main = "",
par.settings = list(axis.line = list(col = 'transparent')))
points(pt2, col = conduct.pal(pt2$conductance_z), cex = 1.5, pch = 20 )
plot(neth_clipped, col = col, col.regions = neth_clipped$group,
edge.col = "transparent", axes = F,
colorkey = list(space = "bottom", height = 0.5, width = 0.7),
main = "",
par.settings = list(axis.line = list(col = 'transparent')))
points(pt3, col = conduct.pal(pt3$conductance_z), cex = 1.5, pch = 20 )
plot(neth_clipped, col = col, col.regions = neth_clipped$group,
edge.col = "transparent", axes = F,
colorkey = list(space = "bottom", height = 0.5, width = 0.7),
main = "",
par.settings = list(axis.line = list(col = 'transparent')))
points(pt4, col = conduct.pal(pt4$conductance_z), cex = 1.5, pch = 20 )
plot(neth_clipped, col = col, col.regions = neth_clipped$group,
edge.col = "transparent", axes = F,
colorkey = list(space = "bottom", height = 0.5, width = 0.7),
main = "",
par.settings = list(axis.line = list(col = 'transparent')))
points(pt5, col = conduct.pal(pt5$conductance_z), cex = 1.5, pch = 20)
plot(neth_clipped, col = col, col.regions = neth_clipped$group,
edge.col = "transparent", axes = F,
colorkey = list(space = "bottom", height = 0.5, width = 0.7),
main = "",
par.settings = list(axis.line = list(col = 'transparent')))
points(pt6, col = conduct.pal(pt6$conductance_z), cex = 1.5, pch = 20)
Skin Conductance by Landcover Types
require(spatialEco)
require(sp)
neth.pts.poly <- point.in.poly(CHIPS_df, neth_clipped)
# Number of points in each landcover group
Conduct_Points_per_group <- tapply(neth.pts.poly@data$conductance_z,
neth.pts.poly@data$group,
FUN=length)
# Mean conductance in each landcover group
Conduct_Mean_per_group <- tapply(neth.pts.poly@data$conductance_z, neth.pts.poly@data$group, FUN=mean)
Conduct_Mean_per_group <- round(Conduct_Mean_per_group, 2)
conductance_groups <- data.frame(t(rbind("Skin Conductance Points (N)"=Conduct_Points_per_group,
"Mean Skin Conductance"=Conduct_Mean_per_group)))
# Number of points in each polygon
Conduct_Points_per_Poly <- tapply(neth.pts.poly@data$conductance_z, neth.pts.poly@data$ITEM2012,
FUN=length)
# Mean conductance in each polygon
Conduct_Mean_per_Poly <- tapply(neth.pts.poly@data$conductance_z, neth.pts.poly@data$ITEM2012, FUN=mean)
Conduct_Mean_per_Poly <- round(Conduct_Mean_per_Poly, 2)
conductance_polys <- data.frame(t(rbind("Skin Conductance Points (N)"=Conduct_Points_per_Poly,
"Mean Skin Conductance"=Conduct_Mean_per_Poly)))library(tableHTML)
# Table: Number of points and means (landcover groups)
conductance_groups[is.na(conductance_groups)] <- 0
conductance_groups_table <- conductance_groups %>%
tableHTML( border = 2,
rownames = TRUE,
headers = c("Sampled Points (N)", "Skin Conductance (M)"),
second_headers = list(c(2, 2), c('Landcover Class', 'Statistic'))) %>%
add_css_second_header(css = list(c('background-color', 'color', 'height'),
c('#C0C0C0', 'black', '50px')),
second_headers = 1:2) %>%
add_css_row(css = list('background-color', '#f2f2f2'), rows = odd(3:6)) %>%
add_css_conditional_column(conditional = ">",
value = 4000,
css = list('background-color', "lightcoral"),
columns = c("Mean Skin Conductance")) %>%
add_css_conditional_column(conditional = "<=",
value = -1000,
css = list('background-color', "lightsteelblue"),
columns = c("Skin Conductance (M)"))
# Table: Number of points and means (all landcover types)
conductance_polys[is.na(conductance_polys)] <- 0
conductance_polys_table <- conductance_polys %>%
tableHTML( border = 2,
rownames = TRUE,
headers = c("Sampled Points (N)", "Skin Conductance (M)"),
second_headers = list(c(1, 3), c('Landcover Class', 'Statistic'))) %>%
add_css_second_header(css = list(c('background-color', 'color', 'height'),
c('#C0C0C0', 'black', '50px')),
second_headers = 1:2) %>%
add_css_row(css = list('background-color', '#f2f2f2'), rows = odd(3:25)) %>%
add_css_conditional_column(conditional = ">",
value = 4000,
css = list('background-color', "lightcoral"),
columns = c("Skin Conductance (M)")) %>%
add_css_conditional_column(conditional = "<=",
value = -1000,
css = list('background-color', "lightsteelblue"),
columns = c("Skin Conductance (M)"))
conductance_groups_table| Landcover Class | Statistic | ||
|---|---|---|---|
| Sampled Points (N) | Skin Conductance (M) | ||
| natural | 33024 | 3250.58 | |
| roads | 46954 | 648.91 | |
| urban | 28333 | 1235.15 | |
conductance_polys_table| Landcover Class | Statistic | ||
|---|---|---|---|
| Sampled Points (N) | Skin Conductance (M) | ||
| Arable land (annual crops) | 6480 | 3134.32 | |
| Construction sites | 0 | 0 | |
| Continuous urban fabric (S.L. : > 80%) | 3995 | 1574.56 | |
| Discontinuous dense urban fabric (S.L. : 50% - 80%) | 4004 | -2737.85 | |
| Discontinuous low density urban fabric (S.L. : 10% - 30%) | 2764 | 5136.59 | |
| Discontinuous medium density urban fabric (S.L. : 30% - 50%) | 1600 | -2537.95 | |
| Discontinuous very low density urban fabric (S.L. : < 10%) | 40 | -0.3 | |
| Fast transit roads and associated land | 0 | 0 | |
| Forests | 9408 | 2292.75 | |
| Green urban areas | 268 | -1840.75 | |
| Herbaceous vegetation associations (natural grassland, moors...) | 2668 | 540.78 | |
| Industrial, commercial, public, military and private units | 15522 | 1851.4 | |
| Isolated structures | 408 | 1944.19 | |
| Land without current use | 1224 | 1316.95 | |
| Mineral extraction and dump sites | 0 | 0 | |
| Open spaces with little or no vegetation (beaches, dunes, bare rocks, glaciers) | 0 | 0 | |
| Other roads and associated land | 43466 | 421.57 | |
| Pastures | 5332 | 867.85 | |
| Permanent crops (vineyards, fruit trees, olive groves) | 340 | 10679.79 | |
| Railways and associated land | 3380 | 3439.65 | |
| Sports and leisure facilities | 7304 | 7481.75 | |
| Water | 108 | 4802.31 | |
| Wetlands | 0 | 0 | |
Space-time Cube
library(OpenStreetMap)
map <- openmap(as.numeric(c(max(pt1$lat), min(pt1$lon))),
as.numeric(c(min(pt1$lat), max(pt1$lon))),
type = "stamen-terrain")
transmap <- openproj(map, projection = proj4string(pt1))
map3d <- function(map, ...){
if(length(map$tiles)!=1){stop("multiple tiles not implemented") }
nx = map$tiles[[1]]$xres
ny = map$tiles[[1]]$yres
xmin = map$tiles[[1]]$bbox$p1[1]
xmax = map$tiles[[1]]$bbox$p2[1]
ymin = map$tiles[[1]]$bbox$p1[2]
ymax = map$tiles[[1]]$bbox$p2[2]
xc = seq(xmin,xmax,len=ny)
yc = seq(ymin,ymax,len=nx)
colours = matrix(map$tiles[[1]]$colorData,ny,nx)
m = matrix(0,ny,nx)
surface3d(xc,yc,m,col=colours, ...) }
library(RColorBrewer)
bp = brewer.pal(11,"RdBu")
library(colourschemes)
cs = rampInterpolate(pt1$conductance_z, rev(bp))
pt1$conduct_pal <- cs(pt1$conductance_z)plot3d(pt1$lon, pt1$lat, pt1$time, xlab="Longitude",
ylab="Latitude", zlab="Time", type = "s",
col = pt1$conduct_pal, size = 2.5)
map3d(transmap)